#delimit;

****;
version 9.0;

clear;
clear matrix;
estimates clear;
eststo clear;
set more off;
set mem 400m;
set seed 365476247;


****************************************;
*STI panel;
cap prog drop runme;
prog def runme;

local hypothesis = `0';
tempfile main bootsave;
use "$da/Crimes-Against-Morality_Main_data-constructed.dta",clear;
keep if dataset==1;

*2. Setting Macros;
*2.1. Setting sample restrictions;

local s3 = "if panelsti_sept==1";
	*in the September STI panel.;

di ;

areg sept_any dd_cl el `s3', cluster(wsid) abs(fe_id);
global mainbeta = _b[dd_cl];
global maint = (_b[dd_cl]-`hypothesis')/_se[dd_cl];
predict epshat, resid;
predict yhat, xb;

*also general the "impose the null hypothesis" yhat and residual;
qui gen temp_y = sept_any - dd_cl*`hypothesis';
areg temp_y el `s3', abs(fe_id);
predict epshat_imposed, resid;
predict yhat_imposed, xb;
qui replace yhat_imposed = yhat_imposed + dd_cl*`hypothesis';

sort wsid el closing;
qui save `main', replace;

*getting a count of then number of clusters in the sample;
qui by wsid: keep if _n==1;
qui summ;
global numclust = r(N);

cap erase `bootsave';
qui postfile bskeep beta_np t_np t_wild
	using `bootsave', replace;

forvalues b=1/$bootreps { ;

qui use `main', replace;
qui by wsid: gen tempbs = uniform();
qui by wsid: gen pos = (tempbs[1] <.5);

qui gen wildresid = epshat_imposed *(2*pos - 1);
qui gen wildy = yhat_imposed + wildresid;
qui areg wildy dd_cl el `s3', cluster(wsid) abs(fe_id);
local bst_wild =(_b[dd_cl]-`hypothesis')/_se[dd_cl];

*nonparametric bootstrap;
bsample, cluster(wsid) idcluster(newsid);
qui areg sept_any dd_cl el `s3', cluster(newsid) abs(newsid);
local bsbeta = _b[dd_cl];
local bst = (_b[dd_cl]-$mainbeta)/_se[dd_cl];

post bskeep (`bsbeta') (`bst') (`bst_wild');
};
*;
qui postclose bskeep;

qui drop _all;
qui set obs 1;
qui gen t_np = $maint;
qui gen t_wild = $maint;
qui append using `bootsave';

qui gen n = .;

foreach stat in t_np t_wild {;
qui summ `stat';
local bign = r(N);
sort `stat';
qui replace n = _n;
qui summ n if abs(`stat'-$maint)< 0.000001;
local myp = r(mean) / `bign';
global pctile_`stat' = 2*min(`myp', (1-`myp'));
};
*;

global mainp = norm($maint);
global pctile_main = 2*min($mainp, (1-$mainp));

local myfmt = "%7.5f";

di;
di "Number BS Reps = $bootreps, Null hypothesis = `hypothesis'";
display "Main beta" _column(13) "main T" _column(22) "Main %le" 
		_column(33) "NP BS %le" _column(44) "wild %le";
di %6.3f $mainbeta _column(13) %6.3f $maint _column(23) 
		`myfmt' $pctile_main _column(34) `myfmt' $pctile_t_np 
		_column(44) `myfmt' $pctile_t_wild ;
		
end;

global bootreps = 999;
runme 0;

